Analyse spatiale et territoriale du logement social
Formation Carthageo-Geoprisme 2022
Introduction
Données RP 2018
Définir le sujet
Soit le sujet : Logements sociaux et qualification des chefs de ménages
Définir les “logements sociaux”
Logements HLM ? Logements SRU ?
Définir la notion de “qualification” ?
Le diplôme le plus élevé ? le nombre d’années d’étude ?
Définir la date
Année 2018 uniquement ? Résultats du RP 2018 (2016-2020) ?
Formuler des questions ou des hypothèses
Qu’elles soient justes ou fausses, les hypothèses permettent de cadrer l’analyse.
Organiser le travail
Sutout dans le cadre d’un groupe !
Ne collecter que les données utiles pour répondre aux questions posées
Afin de ne pas être tenté de partir dans toutes les directions
Archiver soigneusement les programmes et les résultats
Afin de pouvoir reproduire ultérieurement les analyses sur une autre période, un autre territoire
Ne pas attendre d’avoir accumulé tous les résultats pour les commenter
Car l’analyse peut suggérer des erreurs ou ouvrir de nouvelles pistes.
Partir des questions et non pas des outils
Faute de quoi on va trouver des réponses (42 …) sans savoir quelle est la question.
Charger les données statistiques
programme
tab_ind<-readRDS("data/menag2018.RDS")résultat
FALSE COMMUNE ARM IRIS ACHL AEMM
FALSE 1 75056 75101 751010101 C114 2014
FALSE 2 75056 75101 751010101 A11 2008
Préparation de l’analyse
Soit la relation entre logement en HLM (Y) et Diplôme le plus élevé du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés
Le choix de la population de référence est important. Ici on va sélectionner les ménages dont la personne de référence est âgée de 25-39 ans
la sélection ou le regroupement des diplômes est également important car cela va influer sur les résultats du test.
la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage
Sélection des individus et des variables
programme
#table(tab_ind$AGEMEN8)
tab_sel<- tab_ind %>%
filter(AGEMEN8 == "25") %>%
select(DIPLM,HLML, IPONDL) résultats
| DIPLM | HLML | IPONDL |
|---|---|---|
| 18 | 2 | 2.628217 |
| 18 | 2 | 2.628217 |
| 18 | 2 | 2.639089 |
| 18 | 2 | 2.574046 |
Recodage des modalités
On cherche le code des modalités CS1 ezt HLML dans le fichier des métadonnées
meta<-readRDS("data/menag2018_meta.RDS")
metasel <- meta %>% filter(COD_VAR %in% c("DIPLM", "HLML"))
kable(metasel[,c(1,3,4)])| COD_VAR | COD_MOD | LIB_MOD |
|---|---|---|
| DIPLM | 01 | Pas de scolarité ou arrêt avant la fin du primaire |
| DIPLM | 02 | Aucun diplôme et scolarité interrompue à la fin du primaire ou avant la fin du collège |
| DIPLM | 03 | Aucun diplôme et scolarité jusqu’à la fin du collège ou au-delà |
| DIPLM | 11 | CEP (certificat d’études primaires) |
| DIPLM | 12 | BEPC, brevet élémentaire, brevet des collèges, DNB |
| DIPLM | 13 | CAP, BEP ou diplôme de niveau équivalent |
| DIPLM | 14 | Baccalauréat général ou technologique, brevet supérieur, capacité en droit, DAEU, ESEU |
| DIPLM | 15 | Baccalauréat professionnel, brevet professionnel, de technicien ou d’enseignement, diplôme équivalent |
| DIPLM | 16 | BTS, DUT, Deug, Deust, diplôme de la santé ou du social de niveau bac+2, diplôme équivalent |
| DIPLM | 17 | Licence, licence pro, maîtrise, diplôme équivalent de niveau bac+3 ou bac+4 |
| DIPLM | 18 | Master, DEA, DESS, diplôme grande école niveau bac+5, doctorat de santé |
| DIPLM | 19 | Doctorat de recherche (hors santé) |
| DIPLM | YY | Hors résidence principale |
| HLML | 1 | Logement appartenant à un organisme HLM |
| HLML | 2 | Logement n’appartenant pas à un organisme HLM |
| HLML | Y | Hors résidence principale |
s
On recode les modalités des deux variables en regroupant certaines CSP
programme
tab_sel$HLML<-as.factor(tab_sel$HLML)
levels(tab_sel$HLML)<-c("HLM-O","HLM-N",NA)
tab_sel$DIPLM<-as.factor(tab_sel$DIPLM)
levels(tab_sel$DIPLM) <- c("< BAC","< BAC","< BAC","< BAC","< BAC","< BAC",
"BAC","BAC",
"BAC+123","BAC+123","> BAC+3","> BAC+3",NA)
table(tab_sel$DIPLM)FALSE
FALSE < BAC BAC BAC+123 > BAC+3
FALSE 65376 50352 88531 156494
résultats
| DIPLM | HLML | IPONDL |
|---|---|---|
| > BAC+3 | HLM-N | 2.628217 |
| > BAC+3 | HLM-N | 2.628217 |
| > BAC+3 | HLM-N | 2.639089 |
Création du tableau de contingence non pondéré (FAUX)
La solution la plus simple semble être l’instruction table()
programme
tab_cont<-table(tab_sel$HLML,tab_sel$DIPLM)résultats
| < BAC | BAC | BAC+123 | > BAC+3 | Sum | |
|---|---|---|---|---|---|
| HLM-O | 27517 | 16957 | 18316 | 12692 | 75482 |
| HLM-N | 37859 | 33395 | 70215 | 143802 | 285271 |
| Sum | 65376 | 50352 | 88531 | 156494 | 360753 |
Création du tableau de contingence pondéré (JUSTE)
On pondère avec wtd.table() du package questionr.
programme
library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$HLML,tab_sel$DIPLM,
weights = tab_sel$IPONDL)résultats
| < BAC | BAC | BAC+123 | > BAC+3 | Sum | |
|---|---|---|---|---|---|
| HLM-O | 62959 | 38267 | 41036 | 27175 | 169437 |
| HLM-N | 88580 | 77395 | 166707 | 356126 | 688808 |
| Sum | 151539 | 115662 | 207742 | 383302 | 858245 |
Comparaison des niveaux de dépendance automobile
- Tableau non pondéré … légèrement faux !
| < BAC | BAC | BAC+123 | > BAC+3 | Ensemble | |
|---|---|---|---|---|---|
| HLM-O | 42.1 | 33.7 | 20.7 | 8.1 | 20.9 |
| HLM-N | 57.9 | 66.3 | 79.3 | 91.9 | 79.1 |
| Total | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 |
- Tableau pondéré … juste !
| < BAC | BAC | BAC+123 | > BAC+3 | Ensemble | |
|---|---|---|---|---|---|
| HLM-O | 41.5 | 33.1 | 19.8 | 7.1 | 19.7 |
| HLM-N | 58.5 | 66.9 | 80.2 | 92.9 | 80.3 |
| Total | 100.0 | 100.0 | 100.0 | 100.0 | 100.0 |
Visualisation du tableau de contingence
On choisit l’orientation du tableau et on l’affiche avec plot()
mytable<-wtd.table(tab_sel$DIPLM,tab_sel$HLML,weights = tab_sel$IPONDL)
plot(mytable)Visualisation améliorée du tableau de contingence
Tant qu’à faire, on améliore la figure avec des paramètres supplémentaires :
Test du Chi-deux
Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :
mytest<-chisq.test(mytable)
mytest
Pearson's Chi-squared test
data: mytable
X-squared = 97192, df = 3, p-value < 2.2e-16
Visualisation des résidus
Lorsque la relation est significative, on visualise les cases les
plus exceptionnelles avec mosaicplot( …, shade = T)
Localisation aréale (sf)
Le format sf (spatial features)
La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois
- un tableau de données (l’équivalent du fichier .dbf)
- une géométrie (l’équivalent du fichier .shp)
- une projection (l’équivalent du fichier .prj)
Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou
dans d’autres formats standards comme GeoJson, la première tâche
consiste donc à les convertir au formt sf afin de pouvoir les utiliser
facilement dans R. L’importation se fait à l’aide de l’instruction
st_read en indiquant juste le nom du fichier .shp à
charger. Les autres fichiers (.dbf ou .proj) seront lus également et
intégrés dans l’objet qui hérite de la double classe data.frame
et sf.
Etapes de préparation des données
Dans notre exemple, nous allons suivre les étapes suivantes :
- Préparer les données statistiques par IRIS dans un data.frame
- Charger un fonds de carte par IRIS au format sf
- Effectuer une jointure entre les deux fichiers par le code IRIS
- Sauvegarder le résultat
- Agréger les données statistiques et géométriques par commune
- Sauvegarder le résultat.
Préparer les données statistiques
On importe le fichier des individus :
programme
tab_ind<-readRDS("data/menag2018.RDS")résultat
FALSE COMMUNE ARM IRIS ACHL AEMM
FALSE 1 75056 75101 751010101 C114 2014
FALSE 2 75056 75101 751010101 A11 2008
FALSE 3 75056 75101 751010101 A11 2012
Agréger les données
On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :
programme
tab_long<- tab_ind %>%
filter(HLML != "Y")%>%
group_by(IRIS,HLML)%>%
summarise(NB=sum(IPONDL))résultat
| IRIS | HLML | NB |
|---|---|---|
| 751010101 | 1 | 179.08 |
| 751010101 | 2 | 329.51 |
| 751010102 | 1 | 3.16 |
| 751010102 | 2 | 96.72 |
| 751010103 | 1 | 10.97 |
Pivoter le tableau
Puis on fait “pivoter” le tableau pour l’obtenir en format large :
tab_large <- tab_long %>% pivot_wider(id_cols = IRIS,
names_from = HLML,
names_prefix = "HLM_",
values_from = NB,
values_fill = 0)résultat
| IRIS | HLM_1 | HLM_2 |
|---|---|---|
| 751010101 | 179.08 | 329.51 |
| 751010102 | 3.16 | 96.72 |
| 751010103 | 10.97 | 129.56 |
| 751010201 | 148.72 | 1182.91 |
| 751010202 | 16.34 | 923.22 |
Ajouter de nouvelles variables
On ajoute de nouvelles variables telles que le nombre total de ménage et le % de ménages en HLM :
tab<- tab_large %>% mutate(TOT = HLM_1+HLM_2,
HLM_pct = 100*HLM_1/TOT)résultat
| IRIS | HLM_1 | HLM_2 | TOT | HLM_pct |
|---|---|---|---|---|
| 751010101 | 179.08 | 329.51 | 508.59 | 35.21 |
| 751010102 | 3.16 | 96.72 | 99.89 | 3.17 |
| 751010103 | 10.97 | 129.56 | 140.53 | 7.81 |
| 751010201 | 148.72 | 1182.91 | 1331.63 | 11.17 |
| 751010202 | 16.34 | 923.22 | 939.56 | 1.74 |
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par IRIS.
programme
p <- ggplot(tab) + aes (x = HLM_pct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90, 100)) +
scale_x_continuous("% de ménages en HLM") +
scale_y_continuous("Nombre d'IRIS") +
ggtitle(label = "Distribution des logements sociaux dans Paris & PC",
subtitle = "Source : INSEE, RP 2018")résultat
Charger les données géométriques
On importe le fichier des iris du Val-de-Marne qui est au format sf en ne gardant que les colonnes utiles
programme
map_iris <- readRDS("data/map_iris.RDS")
map_iris<-map_iris[,c(4,5,1,2,7)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")résultat
FALSE [1] "sf" "data.frame"
| IRIS | NOM_IRIS | COM | NOM_COM |
|---|---|---|---|
| 751176511 | Ternes 11 | 75117 | Paris 17e Arrondissement |
| 920240403 | Klock | 92024 | Clichy |
Visualisation du fonds iris avec sf
On peut facilement produire une carte vierge des iris du Grand Paris en faisant un plot de la colonne geometry du fichier sf
plot(map_iris$geometry,col="lightyellow")Jointure des données IRIS et du fonds de carte
programme
map_iris_tab<-merge(map_iris,tab,
by.x="IRIS",by.y="IRIS",
all.x=T,all.y=F)résultat
| IRIS | NOM_IRIS | COM | NOM_COM | HLM_1 | HLM_2 | TOT | HLM_pct | geometry |
|---|---|---|---|---|---|---|---|---|
| 751010101 | Saint-Germain l’Auxerrois 1 | 75101 | Paris 1er Arrondissement | 179.08 | 329.51 | 508.59 | 35.21 | MULTIPOLYGON (((651771.6 68… |
| 751010102 | Saint-Germain l’Auxerrois 2 | 75101 | Paris 1er Arrondissement | 3.16 | 96.72 | 99.89 | 3.17 | MULTIPOLYGON (((651668.7 68… |
| 751010103 | Saint-Germain l’Auxerrois 3 | 75101 | Paris 1er Arrondissement | 10.97 | 129.56 | 140.53 | 7.81 | MULTIPOLYGON (((651565.9 68… |
Sauvegarde du fichier par IRIS
On sauvegarde notre fichier au format .RDS de R
saveRDS(map_iris_tab,"data/map_iris_hlm.RDS")Agrégation statistique + géométriques
Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”
Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.
Agrégation des IRIS en communes
L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries
programme
map_com_tab <- map_iris_tab %>%
group_by(COM, NOM_COM) %>%
summarise(HLM_1=sum(HLM_1,na.rm=T),
HLM_2=sum(HLM_2,na.rm=T)) %>%
st_cast("MULTIPOLYGON")
map_com_tab <- map_com_tab %>% mutate(TOT = HLM_1+HLM_2,
HLM_pct = 100*HLM_1/TOT) résultat statistique
| COM | NOM_COM | HLM_1 | HLM_2 | TOT | HLM_pct |
|---|---|---|---|---|---|
| 75101 | Paris 1er Arrondissement | 952.20 | 8321.07 | 9273.27 | 10.27 |
| 75102 | Paris 2e Arrondissement | 562.85 | 11696.77 | 12259.63 | 4.59 |
| 75103 | Paris 3e Arrondissement | 1517.28 | 18267.52 | 19784.80 | 7.67 |
résultat géométrique
Examiner la distribution statistique
On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par Commune.
programme
p <- ggplot(map_com_tab) + aes (x = HLM_pct) +
geom_histogram(breaks = c(0,10,20,30,40,50,
60,70,80,90,100)) +
scale_x_continuous("% de ménages en HLM") +
scale_y_continuous("Nombre de communes") +
ggtitle(label = "Distribution des logements sociaux dans Paris et PC",
subtitle = "Source : INSEE, RP 2018")résultat
Sauvegarde du fichier par commune
On sauvegarde notre fichier au format .RDS de R
saveRDS(map_com_tab,"data/map_com_hlm.RDS")Visualisation (mapsf)
Le package map_sf
Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.
On trouvera la documentation du package mapsf à l’adresse suivante :
Création d’un template cartographique
Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :
map_iris<-readRDS("data/map_iris.RDS")
map_com <-readRDS("data/map_com.RDS")
map_dep <-readRDS("data/map_dep.RDS")
map_iris_hlm<-readRDS("data/map_iris_hlm.RDS")
map_com_hlm<-readRDS("data/map_com_hlm.RDS")tracé d’un fonds de carte vierge
La fonction mf_map() avec le paramètre
type = "base"permet de tracer une carte vide
programme
mf_map(map_iris, type = "base")résultat
Superposition de couches
On peut toutefois ajouter toute une série de paramètres
supplémentaire (col=, border=,
lwd=, …) et superposer plusieurs fonds de carte avec le
paramètre add = TRUE. L’ajout de la fonction
layout permet de rajouter un cadre une légende.
programme
# Trace les Iris avec des paramètres
mf_map(map_iris, type = "base",
col = "lightyellow", border="gray50",lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com, type = "base",
col = NA,border="red",lwd=0.6,
add = TRUE)
# Ajoute les contours des département
mf_map(map_dep, type = "base",
col = NA,border="red",lwd=2,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Paris & Petite Couronne",
credits = "Sources : IGN et INSEE")résultat
Ajout d’un thème
On peut finalement modifier l’ensemble de la carte en lui ajoutant
une instruction mf_theme() qui peut reprendre des styles
existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”,
“candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”,
“barcelona”) mais aussi créer ses propres thèmes
programme
#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris, type = "base",
border="white",
lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com, type = "base",
col = NA, lwd=1,
add = TRUE)
# Ajoute les contours des département
mf_map(map_dep, type = "base",
col = NA, lwd=2,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Paris & Petite Couronne",
credits = "Sources : IGN et INSEE")résultat
Ajout de texte
On peut ajouter une couche de texte avec la fonction
mf_label(). Par exemple, on va ajouter à la carte
précédente le code insee des communes
programme
mf_theme("agolalight")
mf_map(...)
mf_label(map_com,
var="NOM_COM",
cex=0.4,
col="blue",
overlap = FALSE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et Iris du Val de Marne en 2017",
frame = TRUE,
credits = "Sources : IGN et INSEE")résultat
Carte de stock
Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.
Dans le package mapsf, on réalise ce type de carte à
l’aide de la fonction mf_map()en lui donnant le paramètre
type="prop".
On va tenter à titre d’exemple de représenter la distribution du nombre de ménages ordinaires occupant un logement HLM par IRIS :
Carte de stock minimale
Les instructions minimales sont les suivantes :
programme
# Trace les contours des communes
mf_map(x= map_iris,
type = "base")
# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris_hlm,
type ="prop",
var = "HLM_1",
add=TRUE)Mais le résultat est peu satisfaisant car les cercles sont trop
grands. Il faut en pratique toujours effectuer un réglage de ceux-ci
avec l’instruction inches=
résultat
Mais le résultat est peu satisfaisant car les cercles sont trop
grands. Il faut en pratique toujours effectuer un réglage de ceux-ci
avec l’instruction inches=
Carte de stock habillée
programme
mf_theme("agolalight")
mf_map(map_iris, type = "base",
col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="black",lwd=1,add = TRUE)
mf_map(map_iris_hlm, var = "HLM_1",type = "prop",
inches = 0.05, col = "red",leg_pos = "left",
leg_title = "Nombre de ménages", add=TRUE)
mf_layout(title = "Distribution des logements HLM en 2018",
frame = TRUE,
credits = "Sources : IGN et INSEE")résultat
Carte choroplèthe
Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.
La fonction du package mapsf adaptée aux variables
d’intensité est la fonction mf_map()munie du paramètre
type = "choro".
On va prendre l’exemple du nombre de voitures par ménage.
Carte choroplèthe minimale
Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).
programme
# Carte choroplèthe
mf_map(
x = map_iris_hlm,
var = "HLM_pct",
type = "choro")résultat
Carte choroplèthe habillée
On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.
programme
mybreaks = c(0, 10,20,30,40,50,
60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5),
pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris_hlm, var = "HLM_pct",
type = "choro",
breaks = mybreaks,pal = mypal,
border=NA,
col_na = "gray80",leg_title = "% HLM",
leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com, type = "base", col = NA,
border="black",lwd=1,add = TRUE)
mf_layout(title = "% de ménages en HLM au RP 2018",
frame = TRUE,
credits = "Sources : IGN et INSEE")résultat
Carte stock + choroplèthe (1)
On peut combiner les deux modes cartographiques par superposition :
programme
mf_theme("agolalight")
# Choisit les classes
mybreaks = c(0,5,10,20,40,80,100)
# Trace la carte choroplèthe
mf_map(
x = map_iris_hlm,
var = "HLM_pct",
breaks = mybreaks,
# pal=mypal,
type = "choro",
border="white",
col_na = "gray80",
lwd=0.3,
leg_title = "% ménages",
leg_val_rnd = 0,
)
# Ajoute les cercles proportionnels
mf_map(
x =map_iris_hlm,
var = "HLM_1",
type = "prop",
inches = 0.06,
col = "red",
leg_pos = "right",
leg_title = "Nb ménages",
add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com,
type = "base",
col = NA,
border="black",
lwd=1,
add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les ménages ordinaires en HLM 2018",
frame = TRUE,
credits = "Sources : IGN et INSEE")Résultat
Carte stock + choroplèthe (2)
Mais les cercles dissimuent alors les plages de couleur, aussi on
peut utiliser le type prop_choro qui place la variable
choroplèthe à l’intérieur des cercles
programme
mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens","Reds"))
mf_map(map_iris_hlm, type = "base",
col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris_hlm, var = c("TOT", "HLM_pct"),
inches = 0.06, col_na = "grey", pal=mypal,
breaks = mybreaks, nbreaks = 4, lwd = 0.1,
leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
leg_title = c("nb. ménages", "% HLM"),
add = TRUE)
mf_layout(title = "Les ménages ordinaires en HLM 2018",
frame = TRUE, credits = "Sources : IGN et INSEE")résultat
Données spatiales (RPLS, 2020)
Source
Le répertoire des logements locatifs des bailleurs sociaux (RPLS) a pour objectif de dresser l’état global du parc de logements locatifs de ces bailleurs sociaux au 1er janvier d’une année. Il est alimenté par les informations transmises par les bailleurs sociaux. La transmission des informations pour la mise à jour annuelle du répertoire des logements locatifs est obligatoire. Les données sont ensuite géolocalisées à l’adresse et mis à disposition des utilisateurs sur le site du ministère de la transition écologique
Les fichiers sont disponibles en général par régions mais livrés par départements dans le cas de l’Ile de France. Nous allons utilisé ici le fichier du 1er janvier 2020 accessible à l’adresse suivante
https://www.statistiques.developpement-durable.gouv.fr/le-parc-locatif-social-au-1er-janvier-2020-0
Métadonnées
Le fichier de données brutes au format .csv est accompagné d’un document excel précisant le code des variables et la façon dont elles ont été obtenues.
Spatialisation
Le fichier indique pour chaque logement sa localisation précise en terme d’adresse mais aussi d’étage dans un immeuble. A partir de ces données qualitatives, l’INSEE a procédé à un géocodage qui aboutit à la création de deux champs :
- coordonnées de latitude et longitude non projetées
- coordonnées de position en projection Lambert officielle
Selon les analyses on peut utiliser l’une ou l’autre de ces coordonnées. Mais la meilleur solution consiste à créer un fichier de type sf (spatial features) en coordonnées WGS94 qu’on pourra ensuite reprojeter dans le système de son choix.
Stratégie
Avant toute exploitation du fichier il est fortement recommandé d’analyser en détail les métadonnées et de définir une stratégie d’analyse.
- choisir une première zone d’étude de petite taille et localisée de préférence dans un espace que l’on connaît bien.
- choisir des variables intéressantes dont l’on connaît bien la signification et dont on a analysé en détail les métadonnées
- vérifier la qualité des données en regardant notamment le nombre de valeurs manquantes, le dégré de précision, etc.
- sélectionner des données auxiliaires issues d’autres sources que l’on souhaite croiser avec celles du RPLS en s’assurant de leur compatibilité (espace, temps, définition, …)
- Ajouter les coordonnées spatiales et stocker le résultat dans un fichier de type sf comportant les indications de projection.
Importation du fichier
On importe le fichier enregistré au format RDS et on vérifie sa taille avec dim() et sont ype avec class()
don <- readRDS("data/RPLS2020.RDS")
dim(don)[1] 844302 73
Le tableau comporte 844302 lignes (chacune correspondant à un logement) et 73 variables (décrites dans les métadonnées).
Choix de la zone d’étude
On décide de limiter notre analyse dans un premier temps à quatre communes voisines présentant des profils différents. On peut tester leur profil de respect de la loi SRU en 2019 sur l’application suivante : https://www.ecologie.gouv.fr/sru/
- Bonneuil-sur-Marne (94011): large excédent ( > 25% , pas de pénalité)
- Chennevières-sur-Marne (94019) : léger déficit (22.76%, 58 k€)
- Sucy-en-Brie: déficit (94071) : (19.93% , 150k€ )
- Ormesson-sur Marne (94055) : très fort déficit (2.29%, 665 k€)
On relève leur code INSEE afin de pouvoir faciliter l’extraction des données.
Choix des variables
On va se limiter ici à un très petit nombre de variables
variables de localisation
- result_id : code de l’adresse
- result_label : label de l’adresse
- LIBCOM : nom de la commune
- DEPCOM : code de la commune
- latitude : coordonnées latitude
- longitude: coordonnée longitude
- X : coordonnée projetée (EPSG = 2154)
- Y : coordonnée projetée (EPSG = 2154)
variables thématiques
- CONSTRUCT : année de construction
- SURFHAB : surface habitable en m2
- NBPIECE : nombre de pièces
Extraction du fichier
On applique la double sélection des individus et des variables en
nous servant des fonctions filter()et
select()du package dplyr.on aboutit ici à un fichier de
8139 lignes et 11 variables.
sel <- don %>%
filter(DEPCOM %in% c("94011","94019",
"94071","94055")) %>%
select(result_id, result_label,
DEPCOM, LIBCOM,
latitude,longitude,X,Y,
CONSTRUCT,SURFHAB,NBPIECE)
dim(sel)[1] 8139 11
Recodage et typage
Certaines variables doivent être recodées ou changées de type afin de faciliter leur exploitation ultérieure par R.
sel$DEPCOM <- as.character(sel$DEPCOM)
sel$LIBCOM <- as.factor(sel$LIBCOM)
sel$PLG_IRIS <- paste(sel$DEPCOM,sel$PLG_IRIS, sep = "")
sel$SURFHAB <- as.numeric(sel$SURFHAB)Résumé rapide
On analyse rapidement les variables thématiques choisies
| CONSTRUCT | SURFHAB | NBPIECE | |
|---|---|---|---|
| Min. :1902 | Min. : 13.00 | Min. :1.000 | |
| 1st Qu.:1966 | 1st Qu.: 55.00 | 1st Qu.:2.000 | |
| Median :1969 | Median : 68.00 | Median :3.000 | |
| Mean :1977 | Mean : 64.58 | Mean :3.176 | |
| 3rd Qu.:1992 | 3rd Qu.: 77.00 | 3rd Qu.:4.000 | |
| Max. :2019 | Max. :172.00 | Max. :6.000 |
Sauvegarde du fichier
On sauvegarde le fichier obtenu au format .RDS afin de garder le formatage des variables :
saveRDS(sel,"data/sel_logt.RDS")Localisation spatiale
Retour sur sf
Nous revenons sur le package sf (spatial features)
que nous avons déjà rencontré au moment de la création de cartes
thématiques par IRIS ou communes à l’aide du package
mapsf.
Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des logements sociaux au dessus du fonds de carte des IRIS ou communes.
Mais il pourra aussi servir de base à des cartographies
dynamiques permettant de placer les points sur des réseaux de
rue et plus généralement sur des “tuiles” cartographiques permettant
d’effectur des zoom. On utilisera à cet effet d’autres packages comme
leaflet ou sa version simplifiée mapview.
Données ponctuelles
Nous reprenons le fichier de localisation établi au chapitre précédent et nous ne conservons que 6 variables:
logt <- readRDS("data/sel_logt.RDS") %>%
select(adresse=result_id,
X,Y,
date = CONSTRUCT)| adresse | X | Y | date |
|---|---|---|---|
| 94019_0037_00005 | 666779.5 | 6855840 | 1971 |
| 94019_0037_00001 | 666716.7 | 6855829 | 1971 |
| 94019_0037_00001 | 666716.7 | 6855829 | 1971 |
Données IRIS
Nous chargeons par ailleurs le fichier des IRIS en ne gardant que la zone d’étude :
map_iris <- readRDS("data/map_iris.RDS") %>%
filter(INSEE_COM %in% c("94011","94019",
"94071","94055"))| INSEE_COM | NOM_COM | IRIS | CODE_IRIS | NOM_IRIS | TYP_IRIS | DEPT |
|---|---|---|---|---|---|---|
| 94011 | Bonneuil-sur-Marne | 0105 | 940110105 | Saint-Exupéry | H | 94 |
| 94055 | Ormesson-sur-Marne | 0103 | 940550103 | Centre Sud | H | 94 |
| 94071 | Sucy-en-Brie | 0103 | 940710103 | La Cité Verte | H | 94 |
Agrégation par commune
Rappel : on peut agréger les géométries d’un fonds sf. Ici on va créer le fonds de carte des communes.
map_com <- map_iris %>% group_by(INSEE_COM,NOM_COM) %>%
summarise() %>%
st_cast("MULTIPOLYGON")Vérification de la projection
Nous savons que les coordonnées X,Y du fichier logement sont projetées en EPS 2154. Mais quelle est la projection de notre fonds IRIS ? S’agit-il de la même ?
st_crs(map_iris)$proj4string[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
st_crs(2154)$proj4string[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS
Test de superposition
Programme
par(mar=c(0,0,0,0))
#trace les iris
plot(map_iris$geometry,
col="lightyellow", border="gray70",
lwd=0.2)
# trace les communes
plot(map_com$geometry,
col=NA, lwd=1, add=T)
# ajoute les points
points(x=logt$X,
y=logt$Y,
cex=0.2,
col="red",
pch = 16)Résultat
fichier des adresses
Nous allons maintenant établir un fichier de localisation des adresses en nous servant de l’identifiant unique fourni par l’INSEE.
adr <- logt %>% select(adresse,X,Y) %>%
filter(duplicated(adresse) == F) %>%
filter(is.na(X) ==F,is.na(Y)==F)On constate qu’il n’y a que 652 adresses différentes alors que notre fichier fait état de 8139 logements. Une adresse regroupe donc en moyenne plus de 10 logements (habitat collectif).
Transformation en fichier sf
La transformation de notre fichier initial au format sf est facile à
réaliser avec la fonction st_as_sf() du package sf. Mais il
faut prendre garde de bien préciser le système de projection si l’on
veut pouvoir ensuite l’utiliser.
map_adr <- st_as_sf(adr, coords = c("X","Y"))
st_crs(map_adr)<- 2154
str(map_adr)Classes 'sf', 'data.table' and 'data.frame': 612 obs. of 2 variables:
$ adresse : chr "94019_0037_00005" "94019_0037_00001" "94019_0037_00003" "94019_0037_00007" ...
$ geometry:sfc_POINT of length 612; first list element: 'XY' num 666780 6855840
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
..- attr(*, "names")= chr "adresse"
Agrégation des logements
Notre nouveau fichier sf permet désormais d’effectuer des jointures avec le fichier des logements sociaux. A titre d’exemple on peut désormais compter le nombre de logements par adresse et leur ancienneté moyenne.
programme
logt_by_adr <- logt %>%
group_by(adresse) %>%
summarise(nblog = n(),
datemoy = mean(date))résultat
| adresse | nblog | datemoy |
|---|---|---|
| 31 | 2014 | |
| 94011_0017_00001 | 10 | 1970 |
| 94011_0017_00003 | 10 | 1970 |
| 94011_0017_00005 | 10 | 1970 |
| 94011_0019_00002 | 25 | 1992 |
| 94011_0019_00004 | 24 | 1992 |
| 94011_0022_00001 | 20 | 1966 |
| 94011_0022_00002 | 20 | 1966 |
| 94011_0022_00003 | 20 | 1966 |
| 94011_0022_00004 | 20 | 1966 |
Jointure
On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :
map_logt <- inner_join(logt_by_adr,map_adr) %>% st_as_sf()Cartographie avec mapsf
On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :
programme
mf_theme("agolalight")
mybreaks = c(1900, 1950, 1960, 1970, 1980, 1990, 2000, 2010, 2020)
mypal=brewer.pal(n = 8,name = "Spectral")
mf_map(map_iris, type = "base",
col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base",
col = NA,border="black",lwd=1,add = TRUE)
mf_prop_choro( x = map_logt, var = c("nblog", "datemoy"),
inches = 0.08, col_na = "grey", pal=mypal,
breaks = mybreaks, nbreaks = 4, lwd = 0.1,
leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
leg_title = c("nb. logements", "ancienneté"),
add = TRUE)
mf_layout(title = "Les logements sociaux en 2020",
frame = TRUE, credits = "Sources : IGN et RPLS")résultat
Sauvegarde des fichiers cartographiques
On sauvegarde nos différents fichiers cartographiques au format sf relatifs à la zone d’étude.
saveRDS(map_com,"data/sel_map_com.RDS")
saveRDS(map_iris,"data/sel_map_iris.RDS")
saveRDS(map_logt,"data/sel_map_logt.RDS")Cartographie dynamique
Statique ou dynamique ?
- Cartographie statique
- production d’images fixes de qualité
- respect strict des règles de la sémiologie graphique
- choix libre d’une projection adaptée (e.g. EPSG 2154)
- production de documents imprimés à finalité normative ou scientifiques
- Cartographie dynamique
- production d’interfaces consultables dans un navigateur.
- modification possible de l’échelle et de l’arrière-plan
- projection imposée par les “tuiles” (EPSG 4326)
- production de documents interactifs à finalité citoyenne ou exploratoire
Packages R de cartographie dynamique
- leaflet : la référence
- Une librairie javascript non liée à un langage (R, Python, html, …)
- Disponible dans R sous forme de package
- Développement constant
- ggmap : l’empire contre attaque
- des outils cartogaphiques utilisant la syntaxe de tidyverse
- impose désormais un lien avec Google
- tmap : une solution hybride
- permet de passer facilement du mode statique au mode dynamique
- mapview : l’équivalent de mapsf
- mis au point par des développeurs allemands
- facilite l’usage de leaflet
- en progrès constant (mais instable)
Préparation des données
On charge les fichiers au format sf et on les transforme en projection WGS94 (EPSG=4326), condition indispensable pour ajouter des “tuiles” dynamiques lors des zoom.
map_com <- readRDS("data/sel_map_com.RDS") %>%
st_transform(4326)
map_iris <- readRDS("data/sel_map_iris.RDS") %>%
st_transform(4326)
map_logt <- readRDS("data/sel_map_logt.RDS") %>%
st_transform(4326)Carte par défaut
Mapview produit par défaut une carte dynamique du fichier sf.
mapview(map_logt)Superposition de couches
On peut créer des couches et les aditionner avec ‘+’ :
m1 = mapview(map_com, zcol = "NOM_COM")
m2 = mapview(map_logt)
m1+m2Exemple complet
On va essayer de reproduire la carte statique faite avec mapsf
# Carte des communes
map1 <- mapview(map_com, lwd=1, legend= FALSE,
alpha.regions = 0.1)
# Carte des iris
map2 <- mapview(map_iris,lwd = 0.3, label= "NOM_IRIS",
legend= FALSE, alpha.regions = 0)
# Carte des logements
map3 <- mapview(map_logt,
zcol = "datemoy",
at = c(1900,1960, 1970,1980,
1990,2000,2010, 2021),
col.regions = brewer.pal(8, "Spectral"),
cex= "nblog")
map1+map2+map3Annexe : Source des données
Cette annexe présente les sources des différentes données et leurs étapes de transformation préalable aux analyses
Données statistiques sur les logements ordinaires en 2018
Nous partirons des fichiers détail de l’INSEE car, à la différence des tableaux prédéfinis, ils permettent virtuellement toutes les formes de croisement d’indicateurs. Ils sont évidemment beaucoup plus volumineux, mais ce sera justement l’occasion pour les étudiants en data mining d’être confrontés à des problèmes d’optimisation et de big data. On trouve leur description détaillée sur le site de l’INSEE
https://www.insee.fr/fr/statistiques/5542867?sommaire=5395764#consulter
Nous avons opté pour le fichier des individus localisés au canton-ou-ville qui présente une grande polyvalence d’usage puisqu’il permet de reconstituer des tableau agrégés ou l’unité de compte peut-être soit le ménage, soit l’individu selon le critère de pondération adopté.
Etape 1 : téléchargement des données et stockage temporaire
Nous allons télécharger ici le fichier des données pour la région Ile-de-France au format .csv et l’enregistrer dans un dossier spécial tmp qui pourra ulétérieurement être détruit ou déplacé afin de libérer de la place
N.B. Ce programme qui prend quelques minutes sera exécuté une seule
fois. On ajoutera ensuite dans l’en-tête du chunk
eval=FALSE ce qui veut dire que ce bloc de code ne sera
plus executé automatiquement lorsqu’on réalise un knit du document Rmd.
Il sera néanmoins toujours possible de l’executer manuellement en
cliquant sur sa petite flèche verte.
library(curl)
### Téléchargement du fichier INSEE
myurl="https://www.insee.fr/fr/statistiques/fichier/5542867/RP2018_LOGEMTZA_csv.zip"
mydestfile = "tmp/menag2018.zip"
curl_download(url=myurl, destfile=mydestfile)
## Decompression du fichier INSEE
unzip(zipfile = "tmp/menag2018.zip",
exdir = "tmp")
## Examen du contenu
list.files("tmp")Nous constatons que le document zippé contenait en fait deux fichiers différents
Le fichier de données individuelles FD_LOGEMTZA_2018.csv : qui pèse au bas mot 522.3 Mo (1 Giga) et dont nous verrons par la suite qu’il comporte 4.3 millions de lignes et 88 colonnes.
Le fichier de métadonnées varmod_LOGEMT_2017.csv : qui ne pèse que 4.6 Mo et comprend la description précise du label de chacune des modalités de variables.
Etape 2 : Transformation des données au format R
L’importation d’un tableau aussi gros (2.8 millions de lignes et 69 colonnes) donne l’occasion de faire quelques tests de vitesses sur les différents packages capables de lire des fichiers .csv.
Nous allons pour cela utiliser la fonction Sys.time()qui
permet de repérer l’heure au début et à la fin d’une action. Les
résultats dépendront évidemment de la vitesse de l’ordinateur. Il s’agit
ici d’un PC récent de puissance moyenne.
Chargement avec la fonction read.csv
- Avec la fonction
read.csvui fait partie du R-base , le temps de chargement est de 30 secondes. Le tableau résultant est de classe data.frame puisque nous avons utilisé une fonction native de R-base
t1<-Sys.time()
tab<-read.csv("tmp/FD_LOGEMTZA_2018.csv", sep = ";", header =T)
t2<-Sys.time()
paste ("chargement effectué en",t2-t1,"secondes")
dim(tab)
class(tab)Chargement avec la fonction read_csv2
- avec la fonction
read_csv2du package readr, le chargement est effectué en 20 secondes sur le même ordinateur. Le tableau résultant garde la classe data.frame mais est aussi un tibble puisque le package readr fait partie de l’écosystème tibble/tidyverse. Le temps de chrgement est donc divisé par deux.
library(readr)
t1<-Sys.time()
tab<-read_csv2("tmp/FD_LOGEMTZA_2018.csv")
t2<-Sys.time()
paste ("chargement effectué en",t2-t1, "secondes")
dim(tab)
class(tab)Chargement avec la fonction fread
- avec la fonction
freaddu package data.table, le chargement est effectué en 6 secondes sur le même ordinateur.Le tableau résultant conserve la classe data.frame mais possède aussi la classe data.table puisque la fonction fread est issue de ce package. Le temps est divisé par cinq comparativement à la fonction de R-base.
library(data.table)
t1<-Sys.time()
tab<-fread("tmp/FD_LOGEMTZA_2018.csv")
t2<-Sys.time()
paste ("chargement effectué en",t2-t1, "secondes")
dim(tab)
class(tab)On voit donc que le temps de chargement peut différer fortement selon le choix des packages. Il en va ensuite de même pour les traitements d’agrégation des données qui seront plus ou moins rapides selon que l’on utilise les fonctions de R-base applicables à un data.frame, celles du package tidyverse applicables à un tibble ou enfin celles du package data.table applicables à un data.table.
Etape 3 : Sélection des données utiles et sauvegarde au format .Rdata
Nos différentes tableaux peuvent être enregistés au format interne de .R ce qui réduira considérablement leur taille par rapport au fichier texte au format csv qui pèse 522 Mo. Nous allons également limiter la taille du document en ne conservant que les données qui nous intéressent, en l’occurence celles des départements de Paris et Petite Couronne.
Comme ces données bvont nous servir durant tout le projet, elles seront stockées dans le dossier data situé à l’intérieur du projet et non pas dans le dossier tmp qui sera détruit si l’on n’en a plus besoin pour libérer de la place.
- N.B. On ramène l’objet à la classe d’objet unique data.frame pour éviter des conflits possibles entre package. On pourra toujours le retransformer ensuite en data.table ou en tibble.
## Chargement avec fread (+ rapide)
tab<-fread("tmp/FD_LOGEMTZA_2018.csv")
## Suppression de la classe data.table
tab<-as.data.frame(tab)
## Selection des données relatives au Val de Marne
sel<-tab %>% mutate(DEPT=substr(COMMUNE,1,2)) %>% filter(DEPT %in% c("75","92","93","94"))
## Vérification des dimensions du tableau
dim(sel)
## Sauvegarde au format RDS
saveRDS(object = sel,
file = "data/menag2018.RDS")On peut effectuer de façon facultative une sauvegarde au format .csv ce qui évitera des problème d’ouverture du fichier .Rdata pour les personnes ayant des versions anciennes de R. Mais du coup cela engendrera un fichier très volumineux (200 Mo).
## Sauvegarde au format CSV (facultatif)
write.table(x=sel,
file = "data/indiv2017.csv",
sep=";",
dec = ".",
fileEncoding = "UTF-8")Etape 4 : Chargement et sauvegarde des méta-données
Il ne faut surtout pas oublier le fichier des métadonnées qui va permettre de recoder facilement tous les facteurs et de décoder les chiffres correspondant aux classes. On va donc le transformer au format R puis l’enregistrer également dans le dossier data.
# Lecture du fichier de métadonnées
meta<-fread("tmp/varmod_LOGEMT_2018.csv")
# Enregistrement dans le dossier data
saveRDS(object = meta,
file = "data/menag2018_meta.RDS")Données géométriques
Les contours des unités spatiales correspondant aux codes de l’INSEE sont produits par l’IGN et disponibles sur le site géoservice en accès libre :
Etape 1 : récupération du fonds IRIS au format shapefile
On récupère le fichier des IRIS et on le décompresse :
https://geoservices.ign.fr/ressource/178706
list.files("tmp/iris")Etape 2 : Importation et transformation au format sf
La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets ubniques rassemblant à la fois
- un tableau de données (l’équivalent du fichier .dbf)
- une géométrie (l’équivalent du fichier .shp)
- une projection (l’équivalent du fichier .prj)
Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou
dans d’autres formats standards comme GeoJson, la première tâche
consiste donc à les convertir au formt sf afin de pouvoir les utiliser
facilement dans R. L’importation se fait à l’aide de l’instruction
st_read en indiquant juste le nom du fichier .shp à
charger. Les autres fichiers (.dbf ou .proj) seront lus également et
intégrés dans l’objet qui hérite de la double classe data.frame
et sf
library(sf)
map <- st_read("tmp/iris/IRIS_GE.shp")
dim(map)
class(map)
head(map,2)Etape 3 : Extraction des IRIS de la zone d’étude
Le fichier comporte près de 50 000 unités spatiales qui correspondent soit à des communes suffisamment grandes pour être découpées en IRIS, soit à des communes non découpées. On reconnaît ces dernières au fait que leur code IRIS se termine par ‘00000’.
Supposons qu’on veuille extraire le fonds de carte du Val de Marne. On va commencer par créer une variable DEPT en extrayant les dxeux premiers caractères du code communal, puis on va sélectionner le départements correspondant :
map_iris<-map %>% mutate(DEPT = substr(INSEE_COM,1,2)) %>%
filter(DEPT %in% c("75","92","93","94"))
dim(map_iris)
class(map_iris)
head(map_iris,2)Le nouveau tableau ne comporte plus que 2752 unités spatiales et 8 colonnes au lieu de 7 puisqu’ l’on a ajouté une colonne DEPT.
On sauvegarde le résultat dans notre dossier data au format interne de R :
saveRDS(object = map_iris,
file = "data/map_iris.Rdata")Etape 4 : création d’un fonds de carte des communes
Comme nous serons amenés à travailler à plusieurs échelles, nous produisons tout de suite un fonds de carte des communes en utilisant les fonctions d’agrégation du packages sf combinées avec celles de dplyr.
map_com <- map_iris %>% group_by(INSEE_COM) %>%
summarise(NOM_COM = min(NOM_COM)) %>%
st_as_sf()
saveRDS(object = map_com,
file = "data/map_com.Rdata")Etape 5 : création d’un fonds de carte par département
Enfin, on construit un fonds de carte des départements selon la même procédure :
map_dep <- map_iris %>% mutate(DEPT = substr(INSEE_COM,1,2))%>%
group_by(DEPT) %>%
summarise() %>%
st_as_sf()
saveRDS(object = map_dep,
file = "data/map_dep.Rdata")Bilan et nettoyage
Nous avons désormais un dossier data qui comporte :
- Le fichier des logements ordinaires en 2018 et ses métadonnées
- Les fonds de carte par iris, commune et département.
- Le fichier du RPLS pour l’année 2020 et sa documentation
list.files("data") [1] "map_com.RDS" "map_com_hlm.RDS" "map_dep.RDS"
[4] "map_iris.RDS" "map_iris_hlm.RDS" "menag2018.RDS"
[7] "menag2018_meta.RDS" "RPLS2018.RDS" "RPLS2020.RDS"
[10] "sel_logt.RDS" "sel_map_com.RDS" "sel_map_iris.RDS"
[13] "sel_map_logt.RDS"
On peut alors décider de détruire le dossier tmp qui contient des dossiers très volumineux et pas forcément indispensables.